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Silicon vacancy center in 4H-SiC: Electronic structure and spin-photon interfaces 
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Defects in silicon carbide are of intense and increasing interest for quantum-based applications due 
to this material’s properties and technological maturity. We calculate the multi-particle symmetry 
adapted wave functions of the negatively charged silicon vacancy defect in hexagonal silicon carbide 
via use of group theory and density functional theory and find the effects of spin-orbit and spin- 
spin interactions on these states. Although we focused on Vg- in 4H-SiC, because of its unique 
fine structure due to odd number of active electrons, our methods can be easily applied to other 
defect centers of different polytpes, especially to the 6H-SiC. Based on these results we identify the 
mechanism that polarizes the spin under optical drive, obtain the ordering of its dark doublet states, 
point out a path for electric field or strain sensing, and find the theoretical value of its ground-state 
zero field splitting to be 68 MHz, in good agreement with experiment. Moreover, we present two 
distinct protocols of a spin-photon interface based on this defect. Our results pave the way toward 
novel quantum information and quantum metrology applications with silicon carbide. 


Over the last several years, deep-center defects in solids 
have been intensely researched for applications in quan¬ 
tum information [1, 2], quantum sensing and nanoscale 
imaging [3] including bioimaging [4, 5]. Their success and 
popularity stem from their unique properties, combining 
advantages from atomic and solid state systems-most no¬ 
tably long spin coherence times even at room tempera¬ 
ture and integrability into a solid state matrix. The NV 
center in diamond is the most studied defect for quan¬ 
tum technologies, so that its properties, strengths and 
limitations are by now very well understood. Deep defect 
centers in silicon carbide (SiC) have emerged as strong 
contenders due to this material’s significantly lower cost, 
availability of mature microfabrication technologies [6, 7], 
and favorable optical emission wavelengths [8]. 

Some of the stable defects in SiC have the same struc¬ 
ture as the NV center in diamond in terms of symmetry 
and the number of active electrons and, as a result, spin 
and electronic structure. Such defects include the silicon- 
carbon divacancy, which has been investigated over the 
last several years [9-12]. Experiments [8, 13-18] on the 
Si monovacancy (VgJ have shown that this is a distinct 
defect in terms of electronic and spin structure. It fea¬ 
tures a ground state with total spin 3/2 [13, 15], offering 
both quantitative improvements and qualitatively new 
capabilities [16] compared to NV-like defects. To date, 
room temperature spin polarization and coherent con¬ 
trol of Vg- have been implemented via electron spin reso¬ 
nance [17, 19] and optically detected magnetic resonance 
(ODMR) [8, 15, 18, 20]. Unlike the well-studied NV cen¬ 
ter in diamond [21-23], theoretical studies of the Vg- in 
SiC have been mostly limited to finding single-particle 
levels and their energies via density functional theory 
(DET) [24-26]. While such DET calculations are an im¬ 
portant first step, it is of crucial importance to obtain 
the multi-particle electronic structure to understand the 


properties of this defect and take full advantage of the 
novel opportunities it affords. 



FIG. 1. (color online) Vg- in 4H-SiC: (a) Csi.-structure of the 
defect, and the optically-active orbitals of Vg- using DFT: (b) 
u {Ai symmetry), (c) v (Ai), and (d) ex,y (E). Only carbons 
near the Vg- are shown for clarity. 

In this Letter we address this need by calculating the 
multi-particle wave functions of Vg- through a combina¬ 
tion of group theory and DFT. We explicitly find the 
ground states as well as the excited state manifolds, con¬ 
sidering both the orbital and the spin degrees of free¬ 
dom. Furthermore, we investigate the effects of spin- 
orbit and spin-spin interactions. Based on these results 
we (i) explain quantitatively the spin polarization mech¬ 
anism in experiments, (ii) find the zero-field splitting, in 
good agreement with experiment, (iii) present a mecha¬ 
nism that allows this defect to be used for electric field 
or strain sensing, and (iv) propose two spin-photon inter¬ 
face protocols enabled by the rich electronic structure of 
this defect, including the generation of strings of entan¬ 
gled photons and the creation of a Lambda system with 
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TABLE L Negatively charged Si vacancy wave functions for various configurations in the three hole representation. The states 
are classified in terms of orbital electronic configuration, total spin (S) and spin projection along the Ca-axis (ms), overall 
symmetry representation of the state (T) and its decomposition in terms of the orbital and spin symmetries (To (^Ts). ql states 
(not explicitly shown) are defined similarly to states Tg-Tg with the replacement v u. The notation || ...) represents the 
Slater determinant of each component inside the bracket. The bar (no bar) over an orbital indicates spin down (up). 


potential applications in quantum technologies. 

The Cezy symmetry of bulk 4H-SiC is lowered to the 
point group in the presence of Vg-. The local ge¬ 
ometry of Vg- is shown in Fig. 1(a), where the missing 
silicon leaves four dangling bonds (sp^-orbitals) on the 
surrounding carbons. Single electron molecular orbitals 
(MO) can be constructed from symmetry-adapted linear 
combinations of the three equivalent sp^-orbitals (a, b 
and c) from the basal-plane carbons and the sp^-orbital, 


d, belonging to the carbon atom on the Cs-axis that coin¬ 
cides with the crystalline c-axis. Using the standard pro¬ 
jection operator technique [27] and our DFT results as a 
guide [Fig. l(b)-(d)], we obtain the following MOs of the 
defect center: u=au{ci-\-b-\-c)-\-Pud^ v=ay{a-\-b-\-c)-\-Pyd, 
ex=Oix{2c—a—b)^ and ey=ay{a—b), where the coefficients 
are given in [28]. The orbitals, as calculated by DFT, are 
shown in Fig.l. The functions u and v transform as Ai, 
ex and ey transform as the x and y components of the 











E representation respectively and the states are listed in 
order of increasing energy according to our DFT calcu¬ 
lations. 

The electronic configuration of this defect is modeled 
by three holes, a simpler but equivalent picture to that of 
five active electrons. Then, the three-hole lowest energy 
quartet configurations are identified as vCxCy, uexCy, and 
uvcx (or uvcy), respectively, increasing in energy [28]. 
The tensor products of u, and Cx^y states with the to¬ 
tal spin eigenstates comprise our basis set, from which 
we calculate the multi-particle symmetry-adapted states 
compatible with Csu- The odd number of particles here 
results in a much more complicated structure compared 
to NV centers in diamond and divacancies in SiC. Thus, 
we obtain the multi-particle wave functions systemati¬ 
cally by use of the projection operator on the basis states 
for both the orbital and the spin degree of freedom: 

( 1 ) 

R 

where, {R) is the character of operation R in the 
irreducible representation [28], and T is the irreducible 
matrix representation for the R symmetry operator (ten¬ 
sor product of the three-particle orbital and spin oper¬ 
ators [28]). The resulting symmetry adapted states are 
shown in Table I, and are characterized by the total spin 
S', the orbital and spin symmetry, as well as their overall 
symmetry. These classifications are of key importance in 
understanding the nature of these states, their additional 
interactions, as well as the allowed optical or spin-orbit 
assisted transitions and selection rules. The ground state 
manifold has S=3/2 (quartet), while there are nearby 
additional manifolds (each a doublet, 5'=l/2) with some 
having the same orbital composition as the ground state 
and split from each other only due to Coulomb interac¬ 
tions (see Fig. 3 and [28]). 

The states are split and mixed further by spin-orbit 
(SO) and spin-spin interactions. The SO coupling is 

= ( 2 ) 

3 

where tj and Sj are orbital and spin angular momen¬ 
tum operators belonging to the hole. The former is 
defined as {^j)i = eiki\VV{rj)]k[pj\i/2rin?c^ where the 
V{rj) is the local potential, pj is the hole momentum 
operator with coordinate indices The compo¬ 

nents of both t and s transform as the {Ey^ Ex, A 2 ) 
representation and the Hso Hamiltonian itself trans¬ 
forms as Ai. With these symmetry classifications we 
see that the diagonal part of Hso^ will 

only couple states of the same L and S and of or¬ 
bital symmetry E (since A1CE0A20E). Thus, the 
ground states do not split due to this term, while states 
{T;^} and {^q 2 } shift and/or mix within their mani¬ 


folds, as shown in Fig. 2 by A(^={(t)^\\Lf^\\(t)^)/{2^/2) 
and Aq=(0f„J|iy||(/)f„J/(2\/2) respectively (given in 
terms of reduced matrix elements and ^ = {e^^v‘^e}). 
Note that the total orbital angular momentum operator 
is used here, which is equivalent to using Eq. 2 for matrix 
elements between states of the same total S and L [27]. 

The transverse parts of the SO interaction, 
couple states of different total spin 
and orbital character {u^v} to both Cx and Cy at single 
particle level. Hence the ground states will couple to 
{^dii (defined in Table I) via these transverse SO terms. 
This coupling is crucial both in explaining existing 
experiments and in designing future applications. The 
key is to notice that ground states and ql excited 
states with \Sz\=3/2 couple more strongly to excited 
{^dii states compared to the states with \Sz\=^/2. 
In fact using the states of Table I we can show that 
the ratio of the matrix elements is v^. From this we 
identify the dominant intersystem crossing channel that 
constitutes the spin polarization mechanism seen in 
recent experiments at the single-spin level [17] with 
h-site (V 2 ) defects, where optical driving polarizes the 
system into the \Sz\=3/2 states. This mechanism, shown 
in Fig. 3, also successfully predicts the recently seen 
increase in the ODMR photo-luminescence intensity 
with microwave drive [8, 15, 17, 18, 20]. 

We can also consider first-order perturbing corrections 
to the ground state wave functions from the excited dark 
doublet states through spin-orbit coupling (see Fig. 3). 
The different strength of the SO matrix elements (e.g., 
the extra involvement of lj,zSj,z with rris = ±1/2 states 
only) will cause a different degree of admixture of excited 
states to the \Sz\=3/2 and \Sz\=l/2 ground states, which 
in turn will allow an electric field [29], strain and mechan¬ 
ical motion [30-32] to couple ground states with different 
\Sz\ projections. This paves the way toward unexplored 
SiC-based applications in sensing. 

Next we consider the spin-spin interaction between the 
holes. The Hamiltonian is 

Rs = E - 3 (si • nj) (sj ■ nj)}, 

i>j b' 

(3) 

where g is the electron g-factor, /tq is the vacuum 
permeability, and (Ib is the Bohr magneton. The 
spin operator of each hole, the distance to each 
other and its unit vector are 5^, Vij and respec¬ 
tively. The spin-spin splittings of the quartets and 
doublets are shown in Fig. 2 in terms of the split- 
ting parameters defined as Zg.='lo{4>te-^\\h\\<i>'^e^)/ 
7qi=7o(0;!:.I \h 7d=7o(<^f 1 1 / 2 110f )/{6^/T0), 
7q2=7o(0f„el 1-^21 |</>f„e)/(2\/l0) and 7q2=7q2(l “ 1-028C), 
where I 2 is an irregular solid harmonic of second rank, 
i.e. 7o=/uo5^/i|/47r, and 

C=(<(>««ell-^2||<(>f„e)/^q~0> see [28]. Using in these expres- 
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FIG. 2. (color online) Electronic configuration of Vg-, shown 
in terms of the wave functions given in Table 1. The split¬ 
tings are shown explicitly for the SO and spin-spin interac¬ 
tions. The spin quartets are grouped on the left half whereas 
the metastable doublets are on the right. The states with 
subscript q and d denote excited quartet and doublet states, 
respectively. The dashed (green) arrows indicate the mixing 
due to spin-spin interactions. 


sions the calculated bond lengths (i=2.058 A, a=2.055 
A, and 6>o=35.26° from our DFT results, we estimate 
the zero field splitting (ZFS) to be 2|7g|=68 MHz, in 
good agreement with experiments [13, 15, 18, 20]. How¬ 
ever, we found a negative D for the ground state, i.e. 
Hs—DlS^ — S{S + l)/3], causing = ±1/2 to be en¬ 
ergetically higher than the rris = ±3/2 states contrary 
to the some assumptions of D>0 in literature. In the 
limit of perfect tetrahedral {Td) symmetry, our calcula¬ 
tion also leads to a vanishing ZFS (0 MHz) consistent 
with the lack of any ZFS with Vg- centers in 3C-SiC. 

Based on Table I, the rich structure of the various 
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FIG. 3. (color online) Spin polarization channel of Vg- 
through the spin-orbit assisted dominant intersystem cross¬ 
ing ^A 2 (ue‘^)^ ^E(e^)^^A 2 (ve‘^) and all other allowed chan¬ 
nels are shown in dashed lines. Thicker lines of blue 
and green indicate 3x faster transition rate from or to 
rus = ±3/2 states by the transverse component of spin- 
orbit A_l whereas orange represents a channel via the lon¬ 
gitudinal component. Energies of the doublets are or¬ 
dered in terms of the one-particle Goulomb Hamiltonian 
X = W'^hilcj)) and leading many-particle direct integrals, 
i.e. f = / Paa{ri)VR(ri,r 2 )paa{r 2 )dPridPr 2 , of Goulomb re¬ 
pulsion [28]. 


transitions and immunity to all local perturbing electric 
and strain fields (Kramer’s degeneracy) enable the de¬ 
sign of a spin-photon interface for applications in quan¬ 
tum computing and quantum communications. Below we 
propose two such protocols. First consider the ground 
states with |5';2|=3/2, split by a B-field along the Cs 
axis, ± ^g. The excited states of interest are 

which are degen¬ 
erate energy eigenstates after SO has been included 
(Fig. 2); these states have \Sz\=‘i/2^ and since the g- 
factor is the same in ground and excited states [24] they 
split by the same amount as the lower levels. They 
are also the only states which are not coupled to the 
states of ^^2 ql manifold via terms. The 

allowed optical transitions between these sets of states 
are and with right and left circularly 

polarized light respectively. Fig. 4(a). A coherently ex¬ 
cited superposition of the two excited states decays to an 
entangled spin-photon state, |cr±)±|4^“) |cr—). Re¬ 
peating this process produces additional photons, all en¬ 
tangled with the spin and each other, resulting in a mul¬ 
tiphoton Greenberger-Horne-Zeilinger state. Augment- 












































































ing the optical protocol with microwaves can couple the 
ground states and allow the production of a cluster state 
[33], similarly to a proposal for quantum dots [34, 35]. 

Next we consider a B-field perpendicular to the C 3 
axis. This mixes all four ground states, and from these 
we select and along with the excited state 
(all of them given in [28] in terms of the states of Ta¬ 
ble I). Then a A-system can be formed. Fig. 4(b). This 
three-level system can be used in numerous quantum ap¬ 
plications and demonstrations, including coherent popu¬ 
lation trapping [36], optical spin qubit rotations [37, 38] 
and generation of spin-photon entanglement [2, 39] with 
applications in quantum repeaters [40]. 



FIG. 4. (color online) (a) A B-field parallel to the C 3 axis 
enables the creation of two two-level systems with the same 
transition frequency but orthogonal polarizations. Periodic 
coherent pumping followed by spontaneous emission leads to 
strings of entangled photons, (b) A B-field perpendicular to 
the C 3 axis allows for the creation of a Lambda system. 

In summary, we addressed the crucial need of calcu¬ 
lating the multi-particle fine structure of the silicon va¬ 
cancy defect in SiC. Based on the resulting spectrum we 
identified the intersystem crossing channel that polar¬ 
izes the system, found a mechanism to enable quantum 
sensing applications, and proposed two spin-photon in¬ 
terface protocols. Our work opens further opportunities 
in understanding these defects and in implementing novel 
quantum technological applications. 
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SINGLE PARTICLE MOLECULAR ORBITALS 
General remarks 

Unlike the case of NV center in diamond, and other 
similar defects such as the axial divacancy in SiC, al¬ 
though we expect the single charged Si vacancy to still 
have a symmetry, it can still be interpreted as a very 
weakly broken symmetry. This is because all four 
nearest neighbors to the vacancy are carbon atoms with 
very similar distances (differing by ^1%) from the near¬ 
est silicon atom in a perfect crystal, and we do not expect 
this qualitative feature to change much upon removal of 
the silicon atom. We thus anticipate that the form of Ai 
symmetry molecular orbitals will be very close to those 
ofTd: 

^ ( 1 ) 

= a b c — 3d, 

where the normalization has been omitted for brevity. 
Based on the ‘nearly-T^^’ symmetry we also anticipate 
that state v will be near-degenerate with states e (in 
Td they are degenerate). As we show below, our single¬ 
particle molecular orbitals obtained using DET indeed 
confirm these qualitative expectations. Combining re¬ 
sults from DET with analytical calculations we can de¬ 
rive further information, such as the coefficients in the 
molecular orbitals, overlap integrals and on-site Coulomb 
energies. 

First principles calculations 

In order to complement the main group theoretic 
results, density functional theory (DET) was used to 
obtain single particle molecular orbitals (MOs) of the 
charged Si-vacancy center in 4H-SiC. The ordering of 
the defect states is obtained from the calculated Kohn- 
Sham eigenstates around the bandgap of 4H-SiC. The 
spin-polarized calculations were carried out using the 
Quantum-ESPRESSO package [1], within the general¬ 
ized gradient approximation (GGA) [2] of Perdew-Burke- 
Ernzerhof (PBE) [3]. In this work, we report the results 
for the Vg-^ at the h-site in a 6 x 6 x 2 (576-atoms) super¬ 


cell with F-centered 2x2x2 k-point sampling according 
to Monkhort-Pack method. 




(c) ex,y (E-symmetry) 


FIG. 1. (color online) Isosurface plots (5 x 10“^e/a.i^“^) 
for the optically-active minority spin MOs of the negatively 
charged silicon vacancy center Vg- in 4H-SiC: (a) the highest 
occupied orbital, u, (b) the lowest unoccupied orbital, v, and 
(c) the next higher unoccupied orbital, ex,y. 

The large size of the supercell considered here ensures a 
reduction in the defect-defect interactions. This produces 
nearly-fiat defect states that are labeled as u/u (Ai- 
symmetry), v/v (Ai-symmetry)and e/e (E-symmetry). 
Here, the letters with bar overhead represent the minor¬ 
ity spin state, with the excess of three-electrons in the 
majority spin states. The MOs of the defect plotted in 
Figure 1 differ from those obtained with group theoretic 
methods (using symmetry-adapted sp^-orbitals) in that 
they are not restricted to the dangling bonds only. In 
DFT-calculations no such restriction is made and it in¬ 
cludes contributions from other electronic states of the 
crystal as well. Nonetheless, the defect states can be seen 
to be highly localized on the carbon atoms surrounding 
the defect. The majority spin u, is found to be resonant 
with the valence band, while the higher energy defect- 









































states lie in the band gap. This ordering of the defect 
states can be seen in Fig.2. Thus, the DFT-results re¬ 
produce the correct symmetries expected from the group 
theoretic results and provide the ordering of the defect 
states relative to each other. 


Fermi Energy 

V 

—- 7^ 

u 

Maj. Spin 

Min. Spin 


FIG. 2. The energy ordering of the defect-induced majority- 
and minority-spin states. 

In the main text we used group-theoretic approach to 
obtain single-particle MOs from the symmetry-adapted 
linear combinations of the sp^-orbitals belonging to the 
four carbons surrounding the silicon vacancy. However, 
group theory does not yield the relative ordering of states 
with same symmetry, which can be obtained from the 
DFT calculations. In Fig.3, we choose a different isosur¬ 
face (compared to the isosurface plots in Fig.l) to show¬ 
case the bonding- and the anti-bonding characters of the 
Ai symmetry states u and respectively. Thus, DFT 
results can be used to shed light on the relative ordering 
of the MOs qualitatively (bonding vs. anti-bonding) and 
quantitatively (Fig.2). 


Coulomb interaction and overlap integrals 


The Coulomb interaction Hamiltonian can be grouped 
as 14 = '^ij + Si terms of interactions be¬ 

tween different sites (denoted by ij) and on-site (ii) in¬ 
teractions. Therefore, the Schrodinger equation in the 
basis of sp^ dangling bonds [4] takes the form of 


^ '^aa '^ad ^ 

'^ab '^aa '^ab '^ad 

'^ab '^ab '^aa '^ad 

\ '^ab '^ab '^ab '^dd J 


— En 


f 1 Ai Ai A2 ^ 

Ai 1 Ai A2 

Ai Ai 1 A2 

y A2 A2 A2 1 y 


( 2 ) 


in terms of the overlap integrals Ai = / and 

^2 = / between the bonds. For Eq. (2) to 

have non-trivial solutions for each eigenenergy {n = 
u^v^Cx^Cy)^ the following determinant has to be zero. 


'^aa Erjx 
'^ab Eji\\ 
'^ab Eji\\ 

'^ab Eji^2 


'^ab EfiXi 

'^aa -F'n 

"^ab EfiXi 
'^ab EjiX2 


'^ab EfiXi 
'^ab ExiX\ 
'^aa Efi 

'^ab EjiX2 


'^ad EjiX2 
'^ad EjiX2 
'^ad EjiX2 
'^dd Efi 


= 0 


( 3 ) 

Note that the sites a, 6, c are equivalent due to the sym¬ 
metry of the basal plane. The Coulomb interaction be¬ 
tween the sites a and d is roughly equal to that of be¬ 
tween a and 6, i.e. \vad\ = (1 ~ ^^)|'^a6| where d ^ 0, 
since the bond length along the c-axis is only slightly 
distorted from the basal ones as shown by the density 
functional theory calculations. Since all sites a-d have 
carbon atoms, \vdd\ = I'^aal = '^o- Moreover, the off-site 
Coulomb interactions are smaller than the on-site inter¬ 
actions, because of the 1/r dependence of the electro¬ 
static potentials, which can be expressed as \vab\ = e^aal- 
Solutions of Eq. (3), with the realistic assumption 
limd ^ 0, leads to the energies of MOs: 




(b) (Ai-symmetry) 


(a) n (Ai-symmetry) 

FIG. 3. (color online) Isosurface plots (5 x 10“^e/a.i^“^) for 
the optically-active minority spin MOs with Ai-symmetry, 
showing: (a) bonding character of u, and (b) anti-bonding 
character of v. 
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( 4 ) 


up to O(e^). The coefficients k and Ak are given by 


K, —1 H- Ai + e(l — 3 A 2 ), 

Ak = [A^ + 3A2 + e^(6Ai — 6A2 + 4) (5) 

+e (—6A1A2 — 2Ai + 6A2 — 6A2)] ^ , 

in terms of the overlap integrals Ai, A2, and the off-site 
to on-site Coulomb ratio e. In the case of zero over¬ 
lap between the bonds (Ai = A2 = 0), according to 
Eq. (4), the energies Ey and ^ become equal, i.e. 
Ey = Ee^^y = -t;o(l - e) and Ey, = -'i^o(l + 3e), also 
















(9) 


indicating that < Ey. This can be understood as the 
defect’s asymptotic limit to tetrahedral symmetry. 

The true benefit of the above treatment is realized once 
it is used in conjunction with the energies calculated by 
DFT. By using the MO energies obtained by DFT (Fig.2) 
in Eq. (4), we find the previously unknown overlap in¬ 
tegrals, the on-site potential energy and the Coulomb 
ratio of the defect to be Ai = 0.0034, A2 = 0.054, 
vq = 1.177eV, and e = 0.285, respectively. Furthermore, 
we calculate the eigenfunctions satisfying Eq. (3) as, 

u = au{a + 6 + c) -h Pud, 

V = ay{a + 6 + c) -h pyd, 

= ax{ 2 c- a - b), 
ey = ay(a - h), 

with the coefficients obtained as au = 0.523, Pu = 0.423, 
ay = —0.272, py = 0.882, ax = 0.408, and ay = 0.707. 
The coefficients of u and v only slightly differ from the 
readily known coefficients of Td symmetry [5], i.e. au = 
Pu = 0.5, ay = —0.289, and Py = 0.866. Later on, we 
use these coefficients to estimate the zero-field splitting 
of the ground state leading to a remarkable agreement 
with the experimentally measured values. 

Energy order of the doublets 

Due to the many-particle nature of the doublets, 
we cannot obtain the ordering of the states using 
DET, which is an effective single-particle description of 
the system. Therefore, we analyzed the ordering via 
Coulomb Hamiltonian Hc= ^hi ■ Vee{ri,rj) using 

the wave functions of the states given in Table. I 

of the main text. One electron (hole) Coulomb terms are 
included in hi, whereas Vee{ri,rj) = e‘^/{Aneo\ri — rj\) 
is the two-particle Coulomb repulsion potential. Eigen 
values of hi in MOs basis are represented by y and 
can be estimated from DET. Many-particle Coulomb 
integrals are given as j^i=f pu{l)VeePii{ 2 )d^rid^r 2 , 
jim=f Pii{^)VeePmm{‘ 2 )d^rid^r 2 , and 

kim=Jpim{^)VeePim{‘^)d^rid^r 2 . The integrals jfi, 
jim^ and kim are the one-center Coulomb integral, 
two-particle Coulomb repulsion direct and exchange 
integrals, respectively. Charge density is defined as 
Pim{i) = '0/(^)*'^m(^) belonging to the particle in the 
basis of sp^ hybridized dangling bond wave functions 
with l,m={a,h,c,d} and l^m. We obtain the Coulomb 
energies of doublets as, 

Efa =Xe3 + 0.67 j 1 + 2.33ja6 - 0.33A:a6 (7) 

E^e2 =Xve^ + 0-22jaa + l-22ja6 + l-Shjad 

- l.22kab + 0.78fc„d (8) 

=Xve^ + 0-41jaa + l-04jab + ^-^Gjad 


- OMkab - O.lSkad 

=Xve^ + 0.74j2a + O.TOjab + 1.56jad 

-h l.SOkab — 0.78kad (10) 

=Xv^e + O.OOj^Q + 0.61jdd + 0.40ja6 + l-OOjac/ 

- O.SOkab - O.S9kad (11) 

where the relationship ^ j ^ ^ holds and due to 

the nearly Td symmetry of the center charge localization 
on the basal and z-axis carbon atoms are assumed to 
be similar, i.e., Eurthermore, we obtain the 

ground state energy in a similar way: 

E^ = 1.44(ja6 - kab) + 1.56(jad - Kd)- (12) 

Assuming Xv‘^e > Xve‘^ > the ordering of doublets 
becomes Efs, E^^^s, Efg2, E^^^s, and increasing in 

energy. 

THREE-PARTICLE STATES 

Because of the near-degeneracy of state v with states 
Cx and Cy, it is energetically favorable for two electrons to 
occupy the e states instead of paying the energetic cost of 
doubly occupying the only slightly lower in energy state 
V. As a result, in the ground state the occupied states 
are v, Cx and ey, in the three-hole picture. 

Eor the ve^ ground state, the three hole configuration 
space is spanned by 32 (2 (g) 4 (g) 4) basis functions in the 
form of single particle Kronecker products = ({'c} (g) 
{a, p})^{{ex, ey}^{a, ^})(g)({e^, ey}^{a, p}). However, 
consideration of the Pauli exclusion principle discards 8 
of them leaving 24 basis states. 

Moreover, the single particle irreducible matrix repre¬ 
sentations for the cases where the degeneracy lies only in 
the orbital, only in the spin, or in both spaces are simply 
r£;(R)(g)ls, or Ve{R)^Te^/^{R), respec¬ 

tively. Note that the identity matrices are defined as 1^ 
for the orbital and Is for the spin subspace. The explicit 
form of the matrices T{R) are given in Table 1. 

Eor the multi-particle vCxCy ground state con¬ 
figuration, the irreducible matrix representation 
r^^l{R) can be decomposed into its orbital and 
spin components for each particle, i.e. r^^l{R) = 

[(1 ® ® (Ee; ® 0 (Ee; 0 (R)- In 

this form, application of the projection operator [5] on 
each basis function, 

ij 

= (ij/h) E x^^HRrr^TR)fi (13) 

R A 

yields the symmetry adapted basis functions belonging 
to the representation of the ground state. Character 
table of Csu is given in Table H. This gives us a prescrip- 




TABLE L Irreducible matrix representations of E and E 1/2 
for orbital and spin degrees of freedom, respectively, 
given in helicity basis with e = expz27r/3. 


tion for generating all the partners of any basis function 
belonging to a given representation. Further combina¬ 
tions of these symmetry adapted basis functions are then 
formed according to the spin configurations listed in Ta¬ 
ble III to finally obtain all the quartet and doublet wave 
functions of ve^ configuration listed in Table II of the 
main text. The wave functions for the uve excited state 
quartet (q2) are also produced in the same way. 
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TABLE 11. Character table of double group. 3-particle 
coordinate (spin) space belongs to the first (last) three rows. 


SPIN-ORBIT ASSISTED TRANSITIONS 
AMONGST DARK DOUBLET STATES 

We show the spin-orbit coupling matrix elements be¬ 
tween all doublet manifolds dl—d5 in Table IV using 
Eq.l of our manuscript with the symmetry-adapted ba¬ 
sis functions given in Table 1. The spin-orbit coupling 
parameters that are perpendicular and parallel to the 
C 3 axis of the defect are represented by and A^, re- 


DV2 0 £>1/2 (g, £)1/2 

Ts 

S 

rris 


2)3/2 

3/2 

+ 3/2 

+ 1/2 

- 1/2 

- 3/2 

\aaa) 

\aafi) + \afia) + |/3aa) 
/3/3 q!) + \PaP) + a/3/3) 
I/3/3/3) 

£)1/2 

1/2 

+ 1/2 

- 1/2 

\aPa) — \Paa) 

|/3a/3) - |a/3^) 

£)1/2 

1/2 

+ 1/2 

- 1/2 

\aPa) ± \Paa) — 2\aaP) 
/3a/3) + |a/3^) - 2 |^^q;) 


TABLE III. Free space spin configuration of three holes in 
terms of spin up a and down (5 states reduced into irrep. of 
a quartet and two doublets 


spectively. Each element of the matrix is evaluated by 
I \M'so I l^i) where i and j are the wave functions given 
as the row and column headings. We also omit the dark 
doublets, much higher in energy, lying in between the 
excited quartet states ql and q2. These will either tran¬ 
sition to the lowest excited quartet state (ql) or along 
the doublet ladder to the five lower doublet states. The 
key thing to notice here is, as shown in Table IV, all 
doublet states except d5 have spin-orbit assisted allowed 
transitions to the lowest dl doublet. However, d5 dou¬ 
blet can transition to dl through the other doublets in 
between and also has strong transition rate -just like dl- 
into the ground spin rus = ±3/2 states by itself which 
will assist the optical spin polarization process. There¬ 
fore, any other high lying doublet states we omitted in 
this fine structure will follow the general paths shown in 
our manuscript and will not affect the dominant spin- 
polarization channel identified as to be through the dl 
doublet in our manuscript. Note that dl is energetically 
the closest doublet to the ground state (as shown above) 
and it is also the only one connected to the ql quartet 
with a directly allowed spin-orbit assisted transition. 


SPIN-SPIN INTERACTION 

Spherical tensor components 

The spin dipole-dipole operator given in terms of the 
single particle operators Sd = —3 (s* • 

can be expressed as {A±H±C±D±E^±F}, using 
the following spherical tensor components, 

A = —dy/riTs Y2 B = Y2 ± , 

c,D = U' (44 + 44) > 

E,F = 4^^44- (14) 







TABLE IV. Spin-orbit matrix elements amongst the dark doublet states. Spin-orbit parameters of the defect are given by 
and along the C 3 axis and the basal plane of the defect, respectively. 


Orbital parts of A and B terms involving the spherical 
harmonic Y 2 belong to the Ai symmetry, whereas all 
other terms belong to the E symmetry. Since the ground 
state wave functions (Table II of the main text) possess 
an A 2 orbital symmetry, only A and B terms of Eq. (14) 
will cause the zero field spin-splitting of the ground state; 
however, for a q2 excited state with E orbital symmetry 
and corresponding spin symmetries listed in Table II of 
the main text, all terms can contribute to the splitting. 

We first calculate the matrix elements of Sd for each 
wave function listed in Table II of the main text by di¬ 
rect evaluation of its spin components. The remaining 
spatial dependence of the matrix elements can then be 
analyzed through the Wigner-Eckart theorem using the 
spatial components of the spherical tensors listed above. 


Ground state zero-field spin splitting 


In the main text, we report the zero-field spin (ZFS) 
splitting of the ground state in a compact form as 




(15) 


where I 2 is an irregular solid harmonic of second rank, 
i.e. ip = y47i7(2r+I)y’"V+^ and 70 = ^ 1 %/(Att). 

In its open form, it can be written as 


Ig =70 


- ^{ve^ey 
i>j 


yO 

I I 


-'X^y 


), 


(16) 


where can be treated as a pair operator. In 

terms of the direct and exchange integrals [5, 6 ], the ex¬ 
pectation value of any pair operator E is given by, 

{X\\F\\X) = 

i>j (17) 

-{aiaj\f{i,j)\ajai)}, 

where E = /(lj) total pair operator and X 

is the multi-particle antisymmetrized product (Slater de¬ 
terminant) defined as AP [ai(l)a 2 ( 2 ) • • • an{N)]. In the 
case of ground state ZFS, these are /(i, j) = 

X = AP[v{l)e,{2)ey{3)]. 

A quantitative estimate of ZFS splitting can be ob¬ 
tained by switching back to the atomic orbitals, 

7, = 

= Y ['nad{r~d){'^ - 3 cos^ Oad) + Vab{r~b)] > (18) 

where the rjab = 1.443 and rjad = 1.557 are the respective 
weight factors of the expectation value oi'igi- 

nating from total ad and ab pair contributions of MOs 
after evaluating the determinantal multi-particle wave 
functions according to Eq. (17) and using the explicit 
forms of 14, V, and ex,y given in Eq. ( 6 ). This equation 
can also be written in a more familiar form starting from 




the spin dipole-dipole interaction as 



1 — 3 cos^ 0 




^ 5(5 + 1 ) 


, (19) 


where 7 ^ becomes — 5'(5' + l)/3]. So far we as¬ 

sumed all the charge of unpaired electrons is localized on 
the neighboring carbon atoms. However, as previously 
reported [7], only 62.3% of total charge is localized on 
the neighboring carbon atoms, and this yields to a re¬ 
duction of roughly r = (1 — (0.377)^) = 0.858 in 7 ^, i.e. 

Evaluation of Eq. (18) with these weight factors and 
structure parameters calculated via DET, i.e. Vab = 
3.3563A, Vad = 3.3567A, and Oad = 35.259°, as well as 
accounting for the missing charge, results with a ground 
state ZES splitting of 2^g ^ — 68 MHz {D < 0) for an 
h-site Vg- defect in good agreement with the experimen¬ 
tally observed values [ 8 , 9]. 


LAMBDA SYSTEM 


A A-type three-level system can be created by a mag¬ 
netic field transverse to the Cs-axis. Such a field will 
mix states in the ground and excited manifolds in dif¬ 
ferent ways. This is because in the ground state man¬ 
ifold there is a small spin-spin splitting between states 
with \Sz \ =3/2 and \Sz\ = 1/2, whereas the correspond¬ 
ing states in the excited manifold are split by the much 
larger spin-orbit interaction. We assume a weak enough 
magnetic field such that the coupling of the spin states is 
much smaller than the spin-orbit term Ae. This allows 
the eigenstates in the excited manifold to remain in the 
form shown in the main text (without a B-field). There 
are several choices for the composition of the A system. 
Below we present some of these options. In all cases, the 
lower levels are eigenstates of Sx, which in terms of the 
states in Table I of the main text are given by 


g,x — 


^2 ^ 

g,x — 


^3 ^ 

g,x — 


A A 


_^(1+_ 7 - i)*S - - i'D 


^(1 - 


^4 ^ _ 11^2 I _ ^4\ 

g^x — 4 g 4 g^Y8^ ^ ^ ’ 

in descending {Sx) value, 3/2,1/2,-1/2, -3/2 and where, 
for simplicity, we have ignored a small correction from 
the ZES. 

In the first approach for a A system we can select states 
with the same weight of | ttt): e.g., states 
or excited state manifold then we 

select which is defined in the main text and has well- 
defined projection of spin along the z (or Cs) axis, due to 
the suppression of Zeeman mixing originating from the 
large SO interaction. The effect is similar to the selection 
rules in self-assembled quantum dot electron-trion sys¬ 
tems under a Voigt B field. Because of this composition 
of the Lambda system, the polarization of the two tran¬ 
sitions is the same. The frequency however is different, 
and that degree of freedom can be used as the ‘handle’ 
with which to the emitted photon can be manipulated. 

An alternate scheme for a A system is to select as lower 
levels the eigenstates of Sx with eigenvalues -1/2 and 3/2, 
^ and In the excited state manifold the rel- 

evant state is then (^q 2 + ^q 2 )/v ^5 again mixing 
to the states with different spin projection along z has 
been ignored in the excited manifold due to the large SO 
splitting. We note that here the two transitions have the 
same polarization but, unlike the scheme above, differ¬ 
ent dipole moments, originating from the different coef¬ 
ficients of ^g, ^g in the states ^ and 
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